  /***************************************************************************
  *   Copyright (C) 2007 by Alexandre Masselot,,,   *
  *   alex@alex-laptop   *
  *                                                                         *
  *   This program is free software; you can redistribute it and/or modify  *
  *   it under the terms of the GNU Library General Public License as       *
  *   published by the Free Software Foundation; either version 2 of the    *
  *   License, or (at your option) any later version.                       *
  *                                                                         *
  *   This program is distributed in the hope that it will be useful,       *
  *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
  *   GNU General Public License for more details.                          *
  *                                                                         *
  *   You should have received a copy of the GNU Library General Public     *
  *   License along with this program; if not, write to the                 *
  *   Free Software Foundation, Inc.,                                       *
  *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
  ***************************************************************************/
#include <boost/program_options.hpp>


#include "utils/runlcmsmsreadermgf.h"
#include "utils/deltamassbuilder.h"
#include "utils/dbsequencefactoryfasta.h"
#include "utils/peptidemsmsfactory.h"
#include "utils/peptidespectrummatch.h"
#include "utils/peaktheomsmsfactorybasic.h"
#include "utils/aminoacidinstancerbasic.h"
#include "utils/enzymetrypsin.h"
#include "utils/psmnistmspreader.h"
#include "msmscoring/scoringmodelsharedpeakcount.h"
#include "msmscoring/scoringmodeldancik.h"
#include "score/scorerecorderlist.h"
#include "score/scoredbestof.h"
#include "score/scoredefault.h"

using namespace std;
using namespace bol;
using namespace boost::program_options;

ScoringModelSharedPeakCount scoringmodelSPC;
ScoringModelDancik scoringmodelDancik;

double scorePeptide(Peptide *pept, Spectrum<PeakIntensity*>* spectrum);

//fragmentation factory
PeakTheoMsmsFactoryBasic peaktheofactory;
DeltaMass *tolPrec, *tolFrag;


int main(int argc, char **argv){
  RunLcmsms run;
  DBSequenceFactoryFasta fasta;

  options_description desc("Allowed options");
  desc.add_options()
      ("help", "produce help message!")
      ("fasta", value<string>(), "fasta databank for seqeunces")
      ("mgf", value<string>(), "peak list")
      ("msp", value<string>(), "peak list from nist msp file")
      ("tolprec", value<string>(), "tolerance (Da or ppm) on the precursor")
      ("tolfrag", value<string>(), "tolerance (Da or ppm) on the fragment")
      ;
  variables_map vm;
  store(parse_command_line(argc, argv, desc), vm);
  notify(vm);
  if (vm.count("help")) {
    cout << desc << "\n";
    exit(1);
  }

  if(vm.count("mgf")){
    RunLcmsmsReaderMGF reader;
    reader.read(run, vm["mgf"].as<string>());
  }
  PeptideModifier peptmodifierReader;
  if(vm.count("msp")){
    ModificationDictionary dico;
    dico.put(new Modification("Cys_CAM C 57.052"));
    dico.put(new Modification("CAM_core K 43.02504"));
    dico.put(new Modification("Oxidation_HW HW 15.999"));
    dico.put(new Modification("Oxidation_M M 15.999"));

    peptmodifierReader.setModifdico(&dico);
    peptmodifierReader.addModif(dico["Cys_CAM"], true);
    peptmodifierReader.addModif(dico["CAM_core"], true);
    peptmodifierReader.addModif(dico["Oxidation_HW"], true);
    peptmodifierReader.addModif(dico["Oxidation_M"], true);

    PSMNistMSPReader reader;
    reader.setPeptideModifier(&peptmodifierReader);
    reader.open(vm["msp"].as<string>());
    while(PeptideSpectrumMatch *pmatch=reader.nextMatch()){
      run.msmsCollection().push_back(pmatch->spectrum());
    }
    reader.close();
  }
  if(run.msmsCollection().size()==0){
    cerr<<"missing --mgf=... or --msp=.. argument"<<endl;
    exit(1);
  }

  if(vm.count("tolprec")){
    DeltaMassBuilder builder;
    tolPrec=builder.deltamassInstance(vm["tolprec"].as<string>());
  }else{
    cerr<<"missing --tolprec=(xDa|yppm) argument"<<endl;
    exit(1);
  }
  if(vm.count("tolfrag")){
    DeltaMassBuilder builder;
    tolFrag=builder.deltamassInstance(vm["tolfrag"].as<string>());
  }else{
    cerr<<"missing --tolfrag=(xDa|yppm) argument"<<endl;
    exit(1);
  }

  peaktheofactory.setSeries("y,b,b++,y++",2);

  scoringmodelDancik.series_add("b");
  scoringmodelDancik.series_add("b++");
  scoringmodelDancik.series_add("y");
  scoringmodelDancik.series_add("y++");
  scoringmodelDancik.init();
  //camsi DB3
/*  scoringmodelDancik.setScoreValues("y", 1.1, 1.92, -0.45);
  scoringmodelDancik.setScoreValues("b", .75, 1.57, -.33);
  scoringmodelDancik.setScoreValues("b++", .16, .78, -.086);
  scoringmodelDancik.setScoreValues("y++", .30, 1.15, -.014);*/
  //camsi DB2
  scoringmodelDancik.setScoreValues("y", .37, .704, -0.253);
  scoringmodelDancik.setScoreValues("b", .34, .613, -.248);
  scoringmodelDancik.setScoreValues("b++", .149, .315, -.123);
  scoringmodelDancik.setScoreValues("y++", .138, .325, -.109);
  
  run.buildIndex();
  run.msmsCollection().index()->setTolerance(tolPrec);

  AminoAcidInstancerBasic aainstancer;
  aainstancer.setMasses();

  //build a list of peptide out of the passed one and the modif set in the mofifier
  PeptideModifier peptmodifier;
  //  peptmodifier.addModif(new Modification("Cys_CAM C 57"), false);
  //  peptmodifier.addModif(new Modification("Oxidation HMW 16"), true);

  //set the enzyme and a vector to store cleaved peptides
  Enzyme *enzyme=new EnzymeTrypsin();
  vector<Peptide> vcleavedpept;

    //the factory will be set up to produce peptide according to the given modifier (modifier cannot be modified thereafter)
  PeptideMSMSFactory factory;
  factory.setPeptideModifier(&peptmodifier);
  //the vector in which will be stored the msms combination
  vector<Peptide> vmsmspept;

  //instanciate the MetaScoreRecorder
  ScoreRecorderList scoreRecorder;

  //used to keep the best 10 match
  ScoredBestOf bestof(10);

  if(vm.count("fasta")){
    fasta.openFile(vm["fasta"].as<string>());
    DBSequence entry;
    while(fasta.nextSequence(entry)){
      enzyme->cleave(*entry.sequence(), vcleavedpept);
      for(vector<Peptide>::iterator it=vcleavedpept.begin(); it<vcleavedpept.end(); it++){
        vmsmspept.clear();
        factory.createPeptides(*it, vmsmspept);
        for(vector<Peptide>::iterator it2=vmsmspept.begin(); it2<vmsmspept.end(); it2++){
          vector<const SpectrumLcmsmsMassIndexAnchor*>::const_iterator first;
          vector<const SpectrumLcmsmsMassIndexAnchor*>::const_iterator last;
          run.msmsCollection().index()->getAnchors(it2->getMass(), first, last);
          for(vector<const SpectrumLcmsmsMassIndexAnchor*>::const_iterator it=first; it<last; it++){
            peptmodifier.printPeptide(cout, &(*it2));
            cout<<"\t"<<(it2->length()-2)<<"\t"<<(*it)->spectrum()->description()<<"\t";
            double scoreval=scorePeptide(&(*it2), (Spectrum<PeakIntensity*>*)(*it)->spectrum());
            ScoreDefault *score=new ScoreDefault();
            (*score)[0]=scoreval;
            MetaScorePtr scorePtr(score);
            scoreRecorder.add(scorePtr);
          }
        }
      }
    }
  }else if(vm.count("msp")){
    PSMNistMSPReader reader;
    reader.setPeptideModifier(&peptmodifierReader);
    reader.open(vm["msp"].as<string>());
    while(PeptideSpectrumMatch *pmatch=reader.nextMatch()){
      peptmodifierReader.printPeptide(cout, pmatch->peptide());
      cout<<"\t"<<(pmatch->peptide()->length()-2)<<"\t"<<pmatch->spectrum()->description()<<"\t";
      double scoreval=scorePeptide(pmatch->peptide(), (Spectrum<PeakIntensity*>*)pmatch->spectrum());
      ScoreDefault *score=new ScoreDefault();
      (*score)[0]=scoreval;
      PeptideSpectrumMatchPtr psm(pmatch);
      MetaScorePtr scorePtr(score);
      psm->setScore(scorePtr);
      scoreRecorder.add(scorePtr);
      bestof.add(psm);
    }
  }
  cout<<"list of scores"<<endl;
  cout<<scoreRecorder;
  cout<<"list of bestof"<<endl;
  cout<<bestof;
}

/**
 * 
 * @param pept 
 * @param spectrum 
 */
double scorePeptide(Peptide *pept, Spectrum<PeakIntensity*>* spectrum){
  Spectrum<PeakTheoMsms*> theospectrum;
  SpectraAlignment<PeakTheoMsms, PeakIntensity> spalign;

  theospectrum.clear();
  peaktheofactory.computeTheoMasses(*pept, theospectrum, 2, true);
  spalign.spectrumA(&theospectrum);
  spalign.spectrumB(spectrum);
  spalign.align(tolFrag);
  double score=scoringmodelSPC.score(spalign);
  cout<<"\t"<<score;
  score=scoringmodelDancik.score(spalign);
  cout<<"\t"<<score<<"\n";
  return score;
}
